/* Copyright 2007-2009 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra
{
    /// <summary>
    /// A matrix class using dense storage.
    /// </summary>
    public partial class DenseMatrix
    {
        /// <summary>Calculates the L1 norm.</summary>
        /// <returns>The L1 norm of the matrix.</returns>
        public override double L1Norm()
        {
            return SafeNativeMethods.d_matrix_norm('1', Rows, Columns, mData);
        }

        /// <summary>Calculates the Frobenius norm of this matrix.</summary>
        /// <returns>The Frobenius norm of this matrix.</returns>
        public override double FrobeniusNorm()
        {
            return SafeNativeMethods.d_matrix_norm('F', Rows, Columns, mData);
        }

        /// <summary>Calculates the infinity norm of this matrix.</summary>
        /// <returns>The infinity norm of this matrix.</returns>   
        public override double InfinityNorm()
        {
            return SafeNativeMethods.d_matrix_norm('I', Rows, Columns, mData);
        }

        /// <summary>
        /// Adds another matrix to this matrix.
        /// </summary>
        /// <param name="other">The matrix to add to this matrix.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <strong>null</strong>.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Add(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.Rows != Rows || other.Columns != Columns)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            DenseMatrix dense = other as DenseMatrix;
            if (dense == null)
            {
                base.Add(other);
            }
            else
            {
                SafeNativeMethods.d_vector_add(mData.Length, mData, dense.mData, mData);
            }
        }

        /// <summary>
        /// Negates each element of this matrix.
        /// </summary>        
        public override void Negate()
        {
            SafeNativeMethods.d_scale(mData.Length, -1.0, mData);
        }

        /// <summary>
        /// Subtracts another matrix from this matrix.
        /// </summary>
        /// <param name="other">The matrix to subtract.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <strong>null</strong>.</exception>
        /// <exception cref="NotConformableException">If the two matrices don't have the same dimensions.</exception>
        public override void Subtract(Matrix other)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (other.Rows != Rows || other.Columns != Columns)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            DenseMatrix dense = other as DenseMatrix;
            if (dense == null)
            {
                base.Subtract(other);
            }
            else
            {
                SafeNativeMethods.d_vector_subtract(mData.Length, mData, dense.mData, mData);
            }
        }

        /// <summary>
        /// Multiplies each element of this matrix with a scalar.
        /// </summary>
        /// <param name="scalar">The scalar to multiply with.</param>
        public override void Multiply(double scalar)
        {
            SafeNativeMethods.d_scale(mData.Length, scalar, mData);
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null" />.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Columns != other.Rows</strong>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not the this.Rows x other.Columns.</exception>
        public override void Multiply(Matrix other, Matrix result)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Columns != other.Rows)
            {
                throw new NotConformableException("other", Resources.ParameterNotConformable);
            }

            if (result.Rows != Rows || result.Columns != other.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
            {
                Matrix tmp = result.CreateMatrix(result.Rows, result.Columns);
                Multiply(other, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                DenseMatrix denseOther = other as DenseMatrix;
                DenseMatrix denseResult = result as DenseMatrix;
                if (denseOther == null || denseResult == null)
                {
                    base.Multiply(other, result);
                }
                else
                {
                    SafeNativeMethods.d_matrix_multiply(Rows, other.Columns, Columns, mData, denseOther.mData, denseResult.mData);
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix by a vector and returns the result.
        /// </summary>
        /// <param name="rightSide">The vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If rightSide is <strong>null</strong>.</exception>
        /// <exception cref="NotConformableException">If <c>this.Columns != rightSide.Size</c>.</exception>
        public override void Multiply(Vector rightSide, Vector result)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }

            if (Columns != rightSide.Count)
            {
                throw new NotConformableException("rightSide", Resources.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Rows != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (ReferenceEquals(rightSide, result))
            {
                Vector tmp = result.CreateVector(result.Count);
                Multiply(rightSide, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                DenseVector denseOther = rightSide as DenseVector;
                DenseVector denseResult = result as DenseVector;
                if (denseOther == null || denseResult == null)
                {
                    base.Multiply(rightSide, result);
                }
                else
                {
                    SafeNativeMethods.d_matrix_multiply(Rows, 1, Columns, mData, denseOther.Data, denseResult.Data);
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix by a vector and returns the result. The operation is vector * matrix.
        /// </summary>
        /// <param name="leftSide">The vector to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        /// <exception cref="ArgumentNullException">If leftSide is <strong>null</strong>.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Rows != leftSide.Size</strong>.</exception>
        public override void LeftMultiply(Vector leftSide, Vector result)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (Rows != leftSide.Count)
            {
                throw new NotConformableException("leftSide", Resources.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Columns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }
            if (ReferenceEquals(leftSide, result))
            {
                Vector tmp = result.CreateVector(result.Count);
                LeftMultiply(leftSide, tmp);
                tmp.CopyTo(result);
            }
            else
            {
                DenseVector denseOther = leftSide as DenseVector;
                DenseVector denseResult = result as DenseVector;
                if (denseOther == null || denseResult == null)
                {
                    base.LeftMultiply(leftSide, result);
                }
                else
                {
                    SafeNativeMethods.d_matrix_multiply(1, Columns, Rows, denseOther.Data, mData, denseResult.Data);
                }
            }
        }

        /// <summary>
        /// Returns the transpose of this matrix.
        /// </summary>        
        /// <returns>The transpose of this matrix.</returns>
        public override Matrix Transpose()
        {
            DenseMatrix tranposed = new DenseMatrix(Columns, Rows);
            SafeNativeMethods.d_transpose(Rows, Columns, mData, tranposed.mData);
            return tranposed;
        }

        /// <summary>
        /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to pointwise multiply with this one.</param>
        /// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If this matrix and <paramref name="other"/> are not the same size.</exception>
        /// <exception cref="NotConformableException">If this matrix and <paramref name="result"/> are not the same size.</exception>
        public override void PointwiseMultiply(Matrix other, Matrix result)
        {
            if (other == null)
            {
                throw new ArgumentNullException("other");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (Columns != other.Columns || Rows != other.Rows)
            {
                throw new NotConformableException("result");
            }

            if (Columns != result.Columns || Rows != result.Rows)
            {
                throw new NotConformableException("result");
            }

            if (result is DenseMatrix && other is DenseMatrix)
            {
                DenseMatrix denseResult = (DenseMatrix) result;
                DenseMatrix denseOther = (DenseMatrix) other;
                SafeNativeMethods.d_vector_multiply(mData.Length, mData, denseOther.mData, denseResult.mData);
            }
            else
            {
                base.PointwiseMultiply(other, result);
            }
        }

        /// <summary>
        /// An implementation of the BLAS GEMM routine: C:= alpha*op(A) * op(B) + beta*C.
        /// where op is to transpose the matrix or not, and C is this matrix.
        /// </summary>
        /// <param name="alpha">The alpha value.</param>
        /// <param name="beta">The beta balue</param>
        /// <param name="transposeA">if set to <c>true</c>, transpose <paramref name="a"/>.</param>
        /// <param name="transposeB">if set to <c>true</c>, transpose <paramref name="b"/>.</param>
        /// <param name="a">the A matrix.</param>
        /// <param name="b">the B matrix.</param>
        public override void Gemm(double alpha, double beta, bool transposeA, bool transposeB, Matrix a, Matrix b)
        {
            if (a == null)
            {
                throw new ArgumentNullException("a");
            }
            if (b == null)
            {
                throw new ArgumentNullException("b");
            }

            int aRows = a.Rows;
            int bRows = b.Rows;
            int aColumns = a.Columns;
            int bColumns = b.Columns;

            if (transposeA && transposeB)
            {
                if (aRows != bColumns)
                {
                    throw new NotConformableException();
                }
                if (Rows != aColumns && Columns != bRows)
                {
                    throw new NotConformableException();
                }
            }
            else if (transposeA)
            {
                if (aRows != bRows)
                {
                    throw new NotConformableException();
                }
                if (Rows != aColumns && Columns != bColumns)
                {
                    throw new NotConformableException();
                }
            }
            else if (transposeB)
            {
                if (aColumns != bColumns)
                {
                    throw new NotConformableException();
                }
                if (Rows != aRows && Columns != bRows)
                {
                    throw new NotConformableException();
                }
            }
            else
            {
                if (aColumns != bRows)
                {
                    throw new NotConformableException();
                }
                if (Rows != aRows && Columns != bColumns)
                {
                    throw new NotConformableException();
                }
            }

            if (Precision.EqualsWithinDecimalPlaces(0.0, alpha, Constants.DecimalsForComparison) && Precision.EqualsWithinDecimalPlaces(0.0, beta, Constants.DecimalsForComparison))
            {
                Clear();
                return;
            }

            if (a is DenseMatrix && b is DenseMatrix)
            {
                int k = transposeA ? aRows : aColumns;
                int tranA = transposeA ? 1 : 0;
                int tranB = transposeB ? 1 : 0;

                if (ReferenceEquals(this, a) && ReferenceEquals(a, b))
                {
                    DenseMatrix tmpA = new DenseMatrix(a);
                    DenseMatrix tmpB = new DenseMatrix(b);
                    SafeNativeMethods.d_gemm(Rows, Columns, k, tranA, tranB, alpha, tmpA.mData, tmpB.mData, beta, mData);
                }
                else if (ReferenceEquals(this, a))
                {
                    DenseMatrix tmp = new DenseMatrix(a);
                    SafeNativeMethods.d_gemm(Rows, Columns, k, tranA, tranB, alpha, tmp.mData, ((DenseMatrix) b).mData, beta, mData);
                }
                else if (ReferenceEquals(this, b))
                {
                    DenseMatrix tmp = new DenseMatrix(b);
                    SafeNativeMethods.d_gemm(Rows, Columns, k, tranA, tranB, alpha, ((DenseMatrix) a).mData, tmp.mData, beta, mData);
                }
                else
                {
                    SafeNativeMethods.d_gemm(Rows, Columns, k, tranA, tranB, alpha, ((DenseMatrix) a).mData, ((DenseMatrix) b).mData, beta, mData);
                }
            }
            else
            {
                base.Gemm(alpha, beta, transposeA, transposeB, a, b);
            }
        }
    }
}